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Abstract 

By viewing the operations of the Richardson purification algorithm as a dis- 
crete time dynamical process, we propose a method to overcome the instability 
of the algorithm by controlling chaos. We present theoretical analysis and nu- 
merical results on the behavior and performance of the stabilized algorithm. 
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In 1950 Richardson [||] proposed an algorithm to construct the eigenvectors of a hermitian 
matrix when its eigenvalues are given. The basic idea of the method is very simple: start 
from an arbitrary initial vector, and use the eigen equation to eliminate the unwanted 
components in the initial vector. Specifically, to compute the eigenvector \ek) , we apply the 
operations according to the equation 

h) ~n(H-^)i0 (o) >, (i) 

where | <p ') is an arbitrary initial vector which has finite overlap with \e^) and {ei} is the 
set of eigenvalues of the N x N hermitian matrix H. If the operations can be done to 
infinite precision, this procedure gives the desired eigenvector in N — 1 iterations. However, 
it is known that this algorithm is unstable 0, i.e. an initial vector close to the desired 
eigenvector is driven away from it under the operations of the algorithm. To see this, let us 
consider a vector 

|0> = |^>+E<M^>, (2) 

where Si are small compared to 1. In the step of the algorithm when (H — Ej) is applied on 
the resultant vector is proportional to 

|0') = l^)+E^K)> (3) 
i^k 

where 

S> = (4) 
£k - £j 

The above equations show that the fixed point \Ek) of the map (H — Ej) is unstable. The 
exponents characterizing the instability can be as large as In W, where W is the natural band 
width defined by the difference between the largest eigenvalue and the smallest eigenvalue 
devided by the smallest level spacing. If the natural band width W is large, then the 
algorithm is highly unstable. It was proposed that one may make the algorithm work 
better in special cases by replacing Eqn.(|]) by 
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\s k ) ~n( H -^) ni l^ (0) ) (5) 

with rii > 1. However, this does not overcome the instability of the original algorithm. 
Experimentation with this proposal showed that the scope of its applicability is very limited. 
Therefore the algorithm was concluded to have little practical use [0]. 

In this paper, we revisit the instability problem of the Richardson algorithm. By viewing 
the operations of the algorithm as a discrete time dynamical process and using ideas similar 
to controlling chaos f||-[5[|, we are able to devise a method to stabilize the algorithm. The 
key observation is that by dynamically arranging the order of the operations in the product 
rii^£;(H — e,i) ni | (p^), we can use the instability itself to enhance the amplitude of the desired 
eigenvector in the initial vector while suppressing those of the unwanted ones. In the rest of 
the paper, we first present the stabilized algorithm. Following that we present a theoretical 
analysis of the behavior of the stabilized algorithm. Finally we present numerical results 
on the behavior and the the effectiveness of the stabilized algorithm: a randomly generated 
initial vector converges to the desired eigenvector exponentially on the average under the 
operations of the algorithm. 

To map the operations under the algorithm to a discrete time dynamical process, we 
rewrite Eqn.(|l]) as an iteration: 

|0(n+i)) = (H -£<")) \<f>W), (6) 

where the sequence G {ei, 1 < i < N,i ^ k} defines the dynamics. As the discrete time 
n increases, the initial state vector | (jr®) of the dynamical system is transformed to | <n - ) ) by 
the set of matrices {H — £j : < i < N, i ^ k}. These matrices map the N- dimensional unit 
sphere to itself if the vectors are normalized, which is assumed through out our paper. Eqn.'s 
(||) and (HP show that the behavior of these maps around je^) is very complex. There are 
large number of stable as well as unstable directions. The behavior of the system is further 
complicated by rounding errors on a finite precision computer. One of the consequences of a 
finite precision is that the operation in Eqn.(|6]) becomes nonlinear and the system becomes 
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chaotic. To see this, consider the operation of H — e\ on a vector | <fi) on a finite precision 
computer. To be more precise, we define a function 

^, ? (|0» = (H-£O (7) 

where ? is the precision of the computer, i.e. for any x G (— <j, the computer gives 1+x = 1. 
Let | 0) = |£j) + 7 \ej) and 7 be sufficiently small. We have on the one hand 

^ £lS {\4>)) = ^{\e l )) = {e l -e l )\e l ) . (8) 

On the other hand, 

^ iC (| £ ,))+7^ il ,(|^)) = (et-£i){k> + 7%)}, (9) 

where 

7 ' = 7^^- (10) 

£j — £1 

Clearly if the natural band width W is large enough, the term proportional to 7' on the right 
hand side of Eqn.(|9|) may not be negligible. Therefore the function 3? e , )? (| <$)) is nonlinear. 
We emphasize that the nonlinearity of 3? e , lC (| 0)) also makes the order of the operations 
in Eqn.([l|) or Eqn.([|) important. In Fig.|l], we show the typical behavior of the system 
described by Eqn.(||). The lower panel shows the separation C = II I 4>i)~ I 02) || between 
two initially neighboring vectors as a function of the discrete time. The initial separation 
between these two vectors is about 10~ 16 . The upper panel shows the computed Lyapunov 
exponent A, which converges to about 0.23. The matrix used to obtain these results is a 
4096 x 4096 tridiagonal matrix, which is more precisely defined later in the paper. We choose 
k to be 1. The sequence {e^} which determines the dynamics is periodic and is given by 
eW = e P . , where j = n mod (N — 1) and P is a random permutation of the N — 1 indices 
{2, • ■ -, iV}. The irregular behavior of ( in the lower panel of the figure is not due to any 
randomness. It is due to the deterministic chaotic dynamics defined by Eqn.(||) on a finite 
precision computer. The task of stabilizing the algorithm is to determine dynamically 
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from the set of discrete values {£j, i ^ k} so that | (j)^) converges to as n increases. The 
ideas of controlling chaos come into the construction of the sequence {e^}. It is interesting 
to contrast our case with that of the usual problem of controlling chaos 0|| where a set of 
continuous parameters is adjusted in dragging the system back to the desired periodic orbit. 

The iterative procedure to compute the eigenvector using the stabilized algorithm 
is the following. 

A-l Generate a normalized random initial vector | (f)^} and an random array {af^ : 
< af ] < 1}. 

A-2 In the n'th step of the algorithm, the parameter eW is determined by e^ n ' = Ej, 
where j is such that aj = maxj^/ c {a i - ra ' ) }. The normalized vectors | (j)^) and the 
array {c4 } are updated according to the equations 

= (H-eW) | (n) ), (11) 

«™=*U- =j ' (i2) 

where A is the smallest level spacing, 5 < A is a small number comparable to 
the accuracy 5 (defined in Eqn.flT9|)) to which the eigen equation is satisfied. 

A-3 Terminate the iteration if the error parameter a = || (H — £f.) \ <j>^) || 2 }^ is 
smaller than required. 

Let M be the total number of iterations and be the number of occurrences of £j in the 
sequence {e^}. Clearly = 0. Less obviously m; oc M as M — > oo. The latter relation 
comes about for the following reason: A newly eliminated element in the array {a- n ^} is set 
to a small residual value 5/ A in the iteration. It grows exponentially in the later iterations 
and is chosen to be eliminated again by the algorithm. 

To analyze the behavior of the algorithm, we need to obtain the accumulated effect of 
the iterations on the initial vector. For that purpose, we write 



* (0) > = E<f l*>> 



(13) 



where \ei) are the eigenvectors of the matrix H. On a finite precision computer, the eigen 
equation is only satisfied to about the machine precision, i. e. we have 



(H - Si) \si) = 'YjRij \sj) , 



where i?j is the same order as the machine precision q. The resultant vector is 



(14) 



(M) ) oc II(H - a™) ■ £ <f = £ <T ] \si) , 



,(Af) 



(15) 



where 



ci M) =cj o) n(e*-e (B) )=^n(^-^r^ 



(16) 



and the terms higher order in q have been ignored. To show that | <^ M )) indeed converges 
to , we compute the ratio 



< 



.(0) 



.(0) 

-* 

JO) 



c (0) 



n 



(17) 
(18) 



where the parameter 5 is defined by 



<5 = max{|i? M |}. 



(19) 



We remark that: (1) If 



£fe— £i 



<C 1, which usually is the case in practical applications, the 



vanishing ratio r\ k is dominated by 



'. The initial values of {cf^} do not affect 

(0)l 



the convergence of the algorithm. (2) The motivation of introducing the array {a\ } in 



by af\ The 



the stabilized algorithm is to replace the a priori unknown coefficients 
previous comment justifies our so doing. 

To demonstrate the practical effectiveness and to illustrate the numerical behavior of our 
algorithm, we show two sets of numerical results. In the first example, we deal with com- 
puting the eigenvectors of of a tridiagonal matrix after its eigenvalues have been obtained 
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by the usual implicit QL method |2]||. In the second example, we deal with the problem 
of computing the basis vectors in the decomposition of a direct product of irreducible rep- 
resentations of the SU (2) algebra. In this case, the algorithm is especially effective because 
the eigenvalues of the matrix involved, the total angular momentum operator, are known 
exactly due to the algebraic properties of SU(2). 

In the first example, we demonstrate the behavior and the performance of the stabilized 
algorithm by computing the eigenvectors of a tridiagonal matrix of dimension N = 4096. 
The diagonal terms H^i G [—1,1] are randomly generated with a uniform distribution. Its 
next diagonal terms are all 1. The eigenvalues of the matrix are obtained by the standard 
implicit QL procedure. For the particular matrix we used, the smallest level spacing A 
equals 1.67 x 10~ 6 and the natural band width W equals 2.7 x 10 6 . The parameter 5 should 
be the same order as the machine precision, which is usually 10 -13 with double precision. In 
our calculations, we choose it to be 5 = 10~ 10 . We will comment on the effect of the exact 
magnitude of 5 on the behavior of the algorithm later in the paper. 

In Fig.|2] and Fig.^], we show respectively the error parameter a defined by 

v = {^\\\<P {n) )-\e k )\\ 2 } l > (20) 

as a function of n for k = 1, the lowest eigenvalue, and k = 3097, the eigenvalue next to the 
smallest level spacing A. For the eigenvector je^) in the above equation, we use the converged 
eigenvector obtained from our stabilized algorithm according to A-l, A-2, and A-3. For 
comparison, the results for a calculated from a naive implementation of algorithm without 
controlling chaos are also shown as dashed lines in the figures. In the naive implementation 
of the algorithm, the unwanted components of the initial vector are eliminated periodically 
instead of according to the procedure outlined in A-l, A-2, and A-3. The effectiveness 
of the stabilized algorithm is evident: the error parameter a decays exponentially on the 
average as a function of n after an initial transient period. In particular, for k — 1, we 
see that the initial vector converges to the desired eigenvector within machine precision in 
about M = 300 iterations. This is much fewer than the number of iterations N — 1 = 4095 
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required by an implementation of the original algorithm with infinite precision. This shows 
that by dynamically arranging the order of the operations, we can indeed use the instability 
itself to enhance the desired component in the initial vector while suppressing the unwanted 
ones. For k = 3097, the initial vector manages to converge to the desired eigenvector within 
machine precision in about 1.5 x 10 4 ~ 4N iterations. We find that this eigenvector is the 
most difficult to compute because its eigenvalue is next to the smallest level spacing A. This 
is similar to other algorithms for computing eigenvectors. We also notice that a exhibits 
large fluctuations as a function of n, even though its magnitude always remains at around 
10 -10 . Detailed analysis of the operations of the algorithm indicates that these fluctuations 
are due to the instability of the original algorithm. However, we want to emphasize that 
these fluctuations do not affect the computation of the desired eigenvector since we can 
always terminate the iteration once a is smaller than required. 

In Fig||, we show the quantity rf^ defined in Eqn.(|17D as a function of % for k = 3097. 
In the inset, we show clS db function of i. The similarity between the curves in the figure 



show that is indeed dominated by the contribution from 



as asserted above in 



the analysis of the algorithm. The largest magnitude of is very small, about 10" 25 . This 
guarantees the convergence to the desired eigenvector for almost any initial vector | (f>^). 
We also find that the relative magnitude of m; depends on the local level spacing around the 
eigenvalue e^. Our calculations indicate that is large (small) when the local level spacing 
is large (small). 

We have experimented with our algorithm on variety of matrices with different eigen 
spectra structure. Our experience indicates that convergence is maintained even when S/A 
is comparable to 1. We find that the number of iterations M needed for the algorithm to 
compute a converged eigenvector of the low-lying eigenvalues within machine precision does 
not increase much as N increases for any given class of matrices. This means that for these 
eigenvectors our algorithm converges in about aN 2 operation, where a depends on the eigen 
structure of the matrices and is a very slowly varying function of N. In the case of tridiagonal 
matrices with random diagonal elements and for the eigenvector of the lowest eigenvalue, 
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our calculations on matrices of dimensions from 512 to 131072 indicate that a is basically a 
constant around 300. For a typical eigenvector, the algorithm converges in about M = (3N 
iterations, i.e. (3N 3 operations, where j3 is a number comparable to 1. The exactly value of 
/3 depends on the local level spacing. We find that the parameter 5 used in the algorithm 
has effects similar to that of the temperature in an optimization problem by simulated 
annealing. When 5/ A is small, we find that the vector | cf)^) may be trapped around a 
particular vector for a long time before it manages to escape. On the other hand, if 8/ A is 
large, we find that it takes longer for the initial vector to converge and the error parameter 
a fluctuates significantly as a function of n. Our experience also indicates that refreshing 
(replacing af 1 ^ by its reciprocal, for example) the array {af 1 "*} periodically when a has become 
small, 10~ 5 for example, helps to speed up the convergence. We believe this eliminates the 
bias built into the relative values of {a^} according to Eqn.(|P7|) during the iterations. In 
comparison to other traditional methods of computing eigenvectors, the stabilized algorithm 
has a number of advantages. Here we want to compare it to the inverse iteration method, 
the best traditional method of computing eigenvectors when the eigenvalues are known 0. 
It typically takes about N 3 operations for inverse iteration to compute an eigenvector. This 
is comparable to the number of operations needed by our algorithm in the worst cases. 
However, the inverse iteration method needs more computer memory because it needs to 
store the inverse of the matrix H— Sk- If the matrix H is unstructurally sparse, it is clear that 
our algorithm can take advantage of the sparsity but the inverse iteration method cannot 
because the inverse of H — may not be a sparse matrix. 

In the second example, we use the algorithm to construct the basis vectors in the de- 
composition of a direct product of irreducible representations of the algebra SU(2). The 
specific case we studied is to construct the basis vectors of the irreducible representations 
contained in the product ■ • ■ (g> S. This is a generalization of the usual Clebsch- 

N 

Gordan coefficients used to couple two angular momenta together. In our example, what 
we want to do physically is to construct the eigenstates of the total angular momentum 
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operator L of N spinless electrons with orbital angular momentum S. This is very useful 

in block diagonalizing the Hamiltonian of a many electron system [7[]. To be more specific, 

we choose N = 8 and S — \. From the algebraic properties of SU(2), we know that 
21 21 21 

— ® — — = 0©1©---© 56 for identical fermions. Therefore the eigenvalues of 
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the total angular momentum operator L 2 are known exactly: {ei — 1(1 + 1),0 < I < 56}. 
When the z component of the angular momentum is constrained to be L z = 0, the dimension 
of the matrix L 2 is D = 8512. This is much larger than the number of distinct eigenvalues. 
Therefore the eigenvalues are highly degenerate. To construct the basis vectors of the eigen 
subspace Qi of dimension Di, we simply start from D\ randomly generated initial vectors 
{| (ftl^), I = 1, 2, • • •, Di} and apply our stabilized algorithm. It is easy to show that the final 
eigenvectors are linearly independent with probability 1. 

In Fig.([|), we show the error parameter a as a function of n for I = 0. Since the 
eigenvalues are massively degenerate, we have to modify the definition of a to 

- = {^711 l^ } > -f: le^Xe^ | ( ->>|| 2 }*, (21) 
iy i=i 

where {| e z }} is the known basis set of the eigen subspace Qi. The data shown in Fig. ([5]) 
are obtained for I = 0. The dimension of the eigen subspace is D = 31. The dashed 
line shows a obtained from a naive implementation of the original algorithm defined in 
previous discussions. We see that a randomly generated initial vector converges to the 
desired eigenvector within machine precision in about M = 75 iterations. We note that 
fluctuations in a are also present, similar to the first example. We have studied the behavior 
of the algorithm with other values of /, N, and S where the dimension of the matrix L 2 
can be as large as 10 6 . We find that in all cases, the initial vector converges to the desired 
eigenvector within machine precision in a number of iterations fewer than twice of the number 
of distinct eigenvalues of the operator L 2 . We conjecture that this is true for all values of I, 
N, and S. A useful point to note in speeding up the construction of the basis vectors is that 
we may use the same sequence {e^} for the construction of all of the basis vectors. This 
avoids the computation efforts in constructing {e^} for each of the initial vectors. One can 
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show that the basis vectors generated this way are linearly independent with probability 1 
due to the instability. 

We expect the stabilized algorithm to have similar advantages for other compact Lie 
algebras. The reasons are mainly twofold. (1) The algebraic properties of the Lie algebra 
involved usually determine the eigen spectra of the operators involved. The eigenvalues 
are known exactly and do not need to be computed. (2) The dimensions of the irreducible 
representations are usually large. This means that the eigen spectra of the operators involved 
are highly degenerate. Therefore the number of distinct eigenvalues and the natural band 
widths are small. Because of these features, we find the stabilized algorithm to be a very 
powerful tool in fully utilizing the symmetry of the physical system concerned |7| . 

Before concluding, we remark that the algorithm should work for general non-hermitian 
matrices as well. Depending on whether the desired eigenvector is a left or right eigenvector, 
we need to modify the operations of the matrix on the vector | 0( n )) to left or right multipli- 
cations. As a general remark, we find the idea of viewing the operations of an algorithm as 
a dynamical process very intriguing. It seems reasonable to believe that the explorations in 
this direction should be interesting and fruitful. 

In summary, we have discussed a method to stabilize the Richardson purification algo- 
rithm using controlling chaos. We have presented theoretical analysis and numerical results 
to illustrate the behavior and the effectiveness of the stabilized algorithm. We find that 
by dynamically arranging the order of operations in the originally unstable algorithm, we 
are able to use the instability to enhance the desired components of a initial vector while 
suppressing the unwanted ones so that the initial vector converges to the desired eigenvector 
exponentially. 
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FIGURES 

FIG. 1. Generic behavior of the system. Upper panel: Lyapunov exponent. Lower panel: 
separation between two initially neighboring vectors. 

FIG. 2. The error parameter a as a function of the number of iterations for the lowest eigen- 
value. Solid line: with controlling chaos. Dashed line: without controlling chaos. 

FIG. 3. The error parameter a as a function of the number of iterations for the eigenvalue next 
to the smallest level spacing A. Solid line: with controlling chaos. Dashed line: without controlling 
chaos. 

FIG. 4. log 10 rf^ as a function of the eigenvalue index i for k = 3096. Inset shows mi as a 
function of i. 

FIG. 5. The error parameter a as a function of the number of iterations for / = 0. Solid line: 
with controlling chaos. Dashed line: without controlling chaos. 
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